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We consider models of relativistic matter containing sharp interfaces across which the matter 
model changes. These models will be relevant for neutron stars with crusts, phase transitions, or 
for viscous boundaries where the length scale is too short to be modelled smoothly. In particular we 
look at immerical techniques that allow us to evolve stable interfaces, for the interfaces to merge, 
and for strong waves and shocks to interact with the interfaces. We test these techniques for ideal 
hydrodynamics in special and general relativity for simple equations of state, finding that simple 
level set-based methods extend well to relativistic hydrodynamics. 
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^ , Numerical simulations of compact objects such as neutron stars (NSs) are widely used to predict the nonlinear 

■ dynamics of, for example, binary merger or collapse Accurate and detailed simulations are thought necessary for 
\ constructing gravitational wave templates for use in detectors such as LIGO and VIRGO 0, Q . The dynamics and 

■ wavesignal are in large part determined by the structure of the NS, which in turn is determined by the bulk equation 
' of state (EOS) and the detailed microphysics in key, narrow, regions near the surface, crust and internal interfaces. 

. Most current models of NSs used in numerical simulations use single ideal perfect fluids (plus magnetic fields) . When 
qh" extending this model, the important physics may be different in specific regions of the NS (the surface and exterior, 
I , the crust and the core, regions of superconductivity, or local phase transitions). Some effects, such as viscosity, will 
play a role only in localised regions. In all of these situations there will be complex behaviour concentrated in a 
thin interface layer. The question is how these interfaces can be modelled accurately and practically in a numerical 
simulation. 

Whether an interface can be simulated practically depends on the length scales involved, both in the underlying 
physics and in the simulation. In single fluid stars the surface (interface between fluid and vacuum) is genuinely sharp. 
Even in more realistic situations the physical length scales may be extremely short. A more quantitative example 
is given by the Ekman layer in hot NS cores which is relevant for proto NS modelling. In this case the length scale 
' has been estimated by [1] to have the form 5 ^ [ri/{p^)Y/'^ with the viscosity coefficient 77 being highly temperature 
\ dependent, but with 77 ^ IQ^^ — IQ^^ g cm~^ s^^ being reasonable for T > 10® K. Working in cgs units and choosing 

■ 7] ~ 10^^ g cm~^ s~^ as a representative value leads to S ^ lOQ^^^^ cm for p ^ 10^** g cm"-^. 

^\ [ We then need to consider what physical length scales can be resolved in feasible numerical simulations. Focusing on 
. finite difference simulations using current technology, we assume that it is practical to simulate for N ~ 10"'^'^ timesteps 
^ ' and that we wish to cover 10 ms in physical time. Then our finest timestep must be 10~^^ s which means, by the 
• ^ , CFL criteria, that are smallest grid spacing must be ~ 10~^ cm. Standard numerical methods smear interfaces and 
' contact discontinuities with time (0]) with the width of a sharp interface being spread over A^i/(p+i) zones after N 
^ , timesteps, where p is the order of accuracy of the method. Making the optimistic choice of p — A we see that the 
- - ' interface is smeared over 10^ zones with a physical width of ~ 10 cm. Even with these optimistic assumptions the 
numerical simulation will be unable to distinguish numerical and physical effects for NSs with even moderate rotation 
rates. 

This argument suggests that it is impractical and inefficient to model the detailed physics involved in these interface 
regions directly. However, we still need to incorporate these effects as best we can within our numerical simulations. 
To do this we will model these interfaces as being genuinely sharp; that is, they have zero width. The additional 
physics will then be incorporated through the dynamical behaviour of the interfaces, and the boundary conditions 
imposed there. 

In this paper we will use a number of single ideal fluid components separated by sharp material interfaces. We will 
work in full general relativity using ideal hydrodynamics but restrict to 1 + 1 dimensions for simplicity. The model 
and a way to simulate it numerically is outlined in section |TT1 A specific numerical implementation is described in 
section Hm Finally we show a variety of tests illustrating the advantages and limitations of this simple implementation 
in section HVl 

Throughout the paper we use a signature of (— , +, +, +) and work in geometric {G ~ c — Mq = 1) units. 
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II. MODELLING SHARP INTERFACES 
A. The continuum model 

The continuum model that we wish to simulate considers the spacetime as being composed of separate regions 
denoted f2x- On each region the physical matter model (as described by a stress-energy tensor and equations of 
motion) may be different. Each region obeys the Einstein equations and any additional equations of motion required 
for the matter. In the present paper we restrict ourselves to vacuum or to single ideal relativistic perfect fluids with 
the stress-energy tensor 

T^"" ^ pohu^'u" + pg''"' , (2.1) 
where po is the rest-mass density, h is the specific enthalpy defined by 

h=(^l + e+^y (2.2) 

with e the specific internal energy, p the pressure, and u'^ the four-velocity. 

The full spacetime dynamics is described by the union of all the regions fix ■ At and across the boundaries dQx the 
Einstein equations must be satisfied, and from the results of Q we note that this should not introduce discontinuities 
to invariant spacetime quantities, and hence we expect no additional problems from the spacetime. Therefore all that 
remains to describe the model is to give the conditions satisfied by the matter model at the boundaries dflx , and the 
conditions specifying the location of the boundaries. In what follows we shall use a standard 3-1-1 splitting approach, 
so the spacetime is foliated into slices. The spacetime regions fix and their boundaries become spacelike volumes 
within the slice. 

For this paper we consider the simplest possible situation for the boundary and interface conditions. We will 
assume that there is force balance at the boundary and no diffusion through it. We will assume that there is no 
surface tension or other non-ideal effect acting on or within the boundary. This means that the boundary will advect 
with the neighbouring fluid velocities, and the boundary condition will depend on the pressure normal to it. 

B. Numerical model 

As noted above we only need to specify the boundary location and the conditions on the matter model at and 
across the boundary. We will first discuss the conditions applied to the location of the boundary. 

1. Boundary location 

In order to impose a boundary condition at dflx we first need to know where the boundary is. We assume that 
the location is known in the initial slice. We need to find the best method of identifying the location of the boundary 
at any later time. We note that this method must be able to deal with complex changes in topology of any interface, 
particularly in merger simulations. 

The standard technique for dealing with sharp features such as shocks and interfaces in a single perfect fiuid is to 
capture them. In this approach the interface is smeared over a (hopefully small) number of points, with the appropriate 
solution found in the continuum limit. This is incompatible with the approach we are taking here. As well as the 
length scale issues discussed in the introduction, there are serious numerical problems when the model changes at the 
interface and the physics of the mixture are not taken into account, as illustrated in appendix 1X1 

The alternative is to explicitly track or locate the interface. This is a standard problem in numerical relativity, as 
event and apparent horizons must be located in numerical simulations (for a review see Q)- Of particular relevance 
to us are the techniques used for event horizons (see for example Q) where the interface is implicitly tracked and 
complex changes of topology can be dealt with. The techniques used borrow substantially from those in Newtonian 
CFD which are precisely used to deal with the problem we are interested in here (see [H, ITOi] for detailed descriptions) . 

Explicitly, we will use a level set to track the boundary dVlx- That is, we introduce a scalar function (j) defined 
over the entire spacetime. The location of the boundary is implicitly given by those points where (f> vanishes, i.e. by 

dflx = {x : (/)(x) = 0} . (2.3) 

The condition defining the location of the boundary dflx then becomes (after the 3 -f 1 split) an evolution equation 
for the scalar function (j). This evolution equation is arbitrary except at the boundary. The equation can therefore 
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be chosen for modelling simplicity or numerical convenience, provided that it reproduces the correct behaviour at the 
boundary. 

In this paper we consider the boundary as being advected with the fluid, i.e. it will be Lie dragged by the 4- velocity 
u'^. The obvious evolution equation for the level set is 

£u<^ = 0. (2.4) 

However, the evolution of the level set can then lead to complex behaviour for (j) away from the boundary, as noted by 
e-g- [1, 0- This complexity is irrelevant to the physical behaviour. Methods for "renormalizing" the level set (which 
changes the behaviour of cj) away from the boundary) are typically used to avoid problems. In 1 + 1 dimensions, and 
when only a single boundary is present, these problems can be avoided by using a constant vector instead of the fluid 
4- velocity. Obviously this constant vector must be chosen to match the velocity at the boundary, i.e. equation (|2.4p 
becomes 

C^Uj = 0. (2.5) 
2. Boundary condition 

To complete the description of the model we need to describe the conditions imposed at the boundary. There are a 
number of properties that we would like these conditions to satisfy. We need the conditions to be compatible with the 
continuum model. We would like the conditions to approximate, as closely as possible, the additional microphysics. 
We need the boundary conditions to be practical in a numerical simulation; in particular, it must be possible to use 
them without introducing unphysical oscillations due to any discontinuities present. We would like the boundary 
conditions to be as simple as possible. 

The condition that we study in th is p aper is the simplest successful condition proposed in Newtonian CFD; the 
Ghost Fluid condition proposed by It is intended to model the simple situation described above; two ideal 

fluids separated by an interface with no diffusion. It is simple and practical to implement, but is not correct in all 
circumstances; studies such as [l^, [l^l have shown that incorrect shock structures can be computed when strong 
shocks hit the interface between extremely different EOSs. We use the Ghost Fluid boundary condition to illustrate 
the general approach. 

The Ghost Fluid approach considers the fluid on each spacetime region fix separately. To impose boundary 
conditions at dflx the fluid describing the matter within the spacetime region fix is artificially extended through the 
interface into the neighbouring region(s). This is analogous to the imposition of boundary conditions using ghost zones 
or ghost points. The boundary conditions are then imposed by providing values for the artificially extended fiuid. In 
this approach a condition is given on the pressure, the velocity and the entropy, from which all other quantities can 
be recovered (given an EOS). 

The precise conditions imposed at the interface are the continuity of the pressure and the normal velocity, and that 
the flow is isentropic. In principal for a material interface we would expect the pressure and normal component of 
the velocity to be continuous; entropy and the tangential velocity may jump. To follow this as closely as possible the 
pressure and normal velocity component are copied from the fluid in the neighbouring region fly, whilst the entropy 
and tangential velocity are extrapolated. This is illustrated in 1 + 1 dimensions in flgure[T] 

III. IMPLEMENTATION 
A. Formulation 

1. Hydrodynamics - flat space 

In order to best compare to standard relativistic hydrodynamics simulations we use the conserved formulation 
of . In 1 + 1 Minkowski spacetime using Cartesian coordinates and the test-fluid approximation this is written as 



9tq + 9,f(q) =0, (3.1) 
where the conserved variables q and fluxes f are given in terms of the primitive variables introduced in section HIl bv 

q = {D, S^, rf = {pqW, pohW^,, pohW^ -p~ poWf , (3.2a) 
i^{Dv\S,v^+p,{T+p)vn^ , (3.2b) 
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FIG. 1: Illustrating the Ghost Fluid boundary condition in one dimension. The boundary condition is to be applied to fluid 1. 
This is done by creating an artificial "ghost" fluid in the region covered by fluid 2. This ghost fluid obeys the same equations 
of motion and EOS as fluid 1 and is used to implicitly capture the boundary conditions at the interface. The pressure and 
(normal) velocity are copied from fluid 2. The entropy (and, in higher dimensions, the tangential velocity) are extrapolated 
from fluid 1. The same procedure is applied to fluid 2. 



where we have introduced the 3-velocity 



a \u 



(3.3) 



where the shift /?* and lapse a resulting from the 3 + 1 split take the trivial values in this case. We have also introduced 
the Lorentz factor 

W = {1 - v^v'y^^\ (3.4) 
To close the system an equation of state is required. In this paper we use the simple 7-law equation of state 

P(Po,e) = (7 - l)Poe, (3.5) 

where 7 is the ratio of specific heats. Different regions of spacetime may have different equations of state which is 
modelled simply by providing different values of 7. 
The level set equation p.4p becomes 



(3.6) 



This equation is not in conservation law form. Whilst it can be written in conservation law form this would have 
no advantage, as it is only the zeros of the level set function that have meaning. The non-conservative form of 
equation l|3.6p is evolved directly. 

2. Hydrodynamics - spherical symmetry 
In 1 + 1 spherical symmetry we follow e.g. [Tsl. [T^. [iTj in using polar-areal coordinates. The line element is written 



{t, r)dt^ + c? {t, r)dr^ + {dO^ + sin^ 9 dc/)^) 



(3.7) 



The 3-velocity becomes 



The conservation law form is written 



= 7^- (3-8) 



dt (aq) + ^dr (aar^f (q)) = s(q), (3.9) 

where the variables, fluxes and sources are 

q = (D, Sr,Tf =. {poW, pohW%r, PohW^ -p- poWf , (3.10a) 

f = {Dv^-, Srv'- +p,ir+ P) v^ f , (3.10b) 

s = aa(^0,^^iSrV^ + T+p + D) + ^,~^Sr^ , (3.10c) 

where the spacetime quantity m — § (l — is the mass aspect function. The spacetime quantities satisfy the 
equations 

dta — —AnraaSr, (3.11a) 



drU — a 



m 



r 



47rr(r + D)-— , (3.11b) 



arlog(a) = ^Anr {Srv'' + p) + . (3.11c) 

We will evolve a using equation (|3.11ap and solve for the lapse using the slicing condition (|3.11cp . whilst the error in 
the Hamiltonian constraint (j3.11bp will be used to check the accuracy of the simulations. 

To reduce problems near the coordinate singularity at the origin the equations are rewritten as 



dt (aq) + ^9.3 (aar^f (D (q)) + 5. (aaf (q)) = s(q), (3.12) 
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where the revised fluxes and sources are 

f = {Dv\ Srv\ (r + p) v'f , (3.13a) 
f(2) = (0,p,0)^, (3.13b) 

s = aa{{),-—{Srv'- + T+p + D),-—Sr\ . (3.13c) 
The level set equation p.4p becomes 

dt<i) + av''dr(j)^Q. (3.14) 
As in flat space, this non-conservative form is evolved directly. 

B. Numerical methods 

To evolve in time the method of lines is used with a second order Runge-Kutta algorithm. 

1. Hydrodynamics 



To illustrate the utility of the approach and the Ghost Fluid boundary condition we have used standard High 
Resolution Shock Capturing (HRSC) methods to solve the hydrodynamics equations p.ip and p.l2p . The domain is 
split into cells, and the conserved variables updated using 
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Following e.g. [ig] we concentrate on Reconstruction-Evolution methods, where the variables are reconstructed to cell 
boundaries and a Riemann problem (approximately) solved to find the inter-cell flux For reconstruction we use 

slope-limited TVD methods with van Leer's monotized-centred (MC) limiter ([l^), or the PPM method ( [20l. [2l| ) . 
For simplicity we use the HLLE flux formula [l^l in what follows, although we have seen no qualitative difference 
when the Marquina flux formula [1^ is used instead. 



2. Level set 



The evolution of the level set (equation (|3.6p or (|3.14p l cannot use the same numerical methods as the equation is 
not in conservation law form. Instead we shall follow the methods outlined in evolving the level set equation as 
a Hamilton- Jacobi equation. General high-order methods for evolving Hamilton- Jacobi equations can be found, for 
example, in p3 |. 

The general method for Hamilton-Jacobi equations of the form 



that works even when the level set (f) is not differentiable is to solve 

d 



(3.16) 



(3.17) 



'i is the numerical value of at gridpoint Xi 
/_(_ appropriate winded approximations of dx4>- 



H the numerical Hamiltonian, assumed Lipschitz continuous, 



where 
and Z?j 

It is simple to use ENO or WENO methods to find the appropriate operators (explicit expressions for low 
order are given by Q). In the 1-1-1 examples studied here it has been sufficient to use first order accurate upwind 
differencing. We would not expect this to be sufficient in more complex examples or in higher dimensions. 

We have also used the simple Lax-Friedrichs approximation to the Hamiltonian, 



H{u+,u-) = H 



(3.18) 



where a is the maximum absolute value of H over the entire domain. In higher dimensions it may be costly to compute 
a over the whole domain, so a local Lax-Friedrichs method (where the maximum value of a over the numerical stencil 
is used) can be used instead. 



3. Ghost Fluid boundary condition 

Figure [1] has already been given to illustrate how the Ghost Fluid boundary condition is implemented. To be 
explicit, we consider a numerical grid {x^} on which the level set (pi changes sign between points /,/ + 1. We apply 
the Ghost Fluid boundary condition to fluid 1, which covers the domain to the left of the interface, using where 
necessary information from fluid 2, which covers the domain to the right of the interface. We will apply the boundary 
condition to at least p + \ "ghost" points extending from 2:7+1 to xj+p+i, where p is the stencil size of the numerical 
method. The reason for the use of p -I- 1 ghost points is that at least one additional point is needed if the interface 
should move one point to the right during the evolution step. 

The normal component of the velocity and the pressure of fluid 1, w^^-* and p*-^-* respectively, are copied from fluid 2. 
The entropy is extrapolated from fluid 1. In the Ghost Fluid method zero order extrapolation normal to the interface 
is used. In the 1-1-1 dimensional cases here, the extrapolation is a simple copy. Explicitly we have 

j = /+l,...,/+P+l, (3.19a) 

j = /+l,...,/+P+l, (3.19b) 

j = /+l,...,H-p+l. (3.19c) 
All other hydrodynamical quantities are found from the pressure, entropy, and the equation of state. 
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4- The evolution step 

In simple cases such as flat space or where free evolution of spacetime quantities is used, the evolution step is now 
simply a matter of updating both the fluid and level set. However, should key spacetime quantities (such as o for 
the case given in section ITlI A 2p be updated using non-local equations (for example, the constraint (I3.11b[) could be 
used to give a) then the order of the update becomes important. For this reason we ensure that we update the level 
set using the methods of section IIII B 21 before the hydrodynamical quantities can be updated using the methods of 
section IIITbH 

We represent a complete evolution step, assuming a first order time-stepping algorithm, as: 

1. Given data {q"} for the hydrodynamical variables and {0"} for the level set and timestep n. 

2. Find interface locations {a;" : (/) = O} , j = 0, . . . , iV", where N"^ — 1 is the number of interfaces at time t". For 
simplicity we denote the boundaries of the computational domain by and x^„ respectively. 

3. If there is only one genuine interface approximate the velocity v there. 

4. Compute the updated level set (t>"''^^ , either using the interface velocity approximated in the previous step, or 
using the velocity of the fluid at each point. 

5. Compute the new interface locations {x^^^ : = 0} , j = 0, . . . , iV"+^. 

6. For each domain = [x^, j = 0, . . . , iV" - 1: 

(a) Check the domain still exists at the next time step, f^"'*'"'^ 7^ {0}- If it is, go to the next domain. 

(b) Apply the Ghost Fluid boundary condition to both boundaries, extending the domain by a sufficient number 
of points, as given in section UlI B 31 

(c) Compute the inter-cell fluxes on this extended domain. 

(d) Restrict the updated variables to the domain f^"^^ — [x^^^^ ^x^^^l]. 

(e) Compute primitive variables and apply physical or symmetry boundary conditions as required. 
Extending this to higher order in time is trivial with method of lines type methods. 

IV. TESTS 

The tests used here are intended to see whether the Ghost Fluid method works in relativistic situations as it does 
in Newtonian theory. This will indicate whether this is a viable approach for numerical simulations of more complex 
models of NS interactions such as mergers. The simplest test is a check that the method would keep a single non- 
interacting material interface stable. This is provided in appendix [X] as part of an illustration of why capturing an 
interface is not straightforward numerically, as discussed in section IIIB II 

We start by looking at the interaction between strong nonlinear features and interfaces in flat space, using the 
model of section llll A 11 These ensure that no spurious oscillations are introduced and see how accurately the method 
captures transient and long-lived smooth features. We then move on to full GR tests, where we can assess the accuracy 
and effectiveness with simple spherically symmetric NS models, both stable and nonlinearly perturbed. 

A. Flat space tests 

1. Shock-interface interaction 

Two tests are considered. In both the domain is x G [0, 1]. In the flrst a stable material interface is set up in 
the middle of the domain across which the equation of state changes due to a change in 7. A mild shock is injected 
from the left edge of the domain. By time t = 1 the shock has interacted with the interface, being reflected and 
transmitted. In this test, which is a modified relativistic analogue of test B of [ll|, the precise initial data is 

po = 1.3346, v= 0.1837, p= 1.5, 7= 1.4 x < 0.05 

po = 1, = 0, P= 1, 7= 1.4 0.05 < a; < 0.5 (4.1) 

Pq ~ 1, V = 0, p — 1, 7= 1.67 a; > 0.5. 
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FIG. 2: Results from an isolated shock hitting a contact discontinuity, ('test B' from The transmitted shock is very 

obvious in the velocity and pressure profiles, and convergence of the solution at this point can be seen. The refiected rarefaction 
wave is considerably smaller, but can still be seen. The contact discontinuity, shown by the dashed line, is captured correctly. 
TVD reconstruction using the MC limiter is used, and only 100 points are shown in each plot. 



Results from this test are shown in figure [21 For simphcity we use TVD reconstruction with the HLLE Riemann 
solver. The transmitted shock and reflected rarefaction are clearly seen, with the results converging towards the 
exact solution. The material interface has propagated to the right due to the interaction, and this is also captured 
accurately. No spurious oscillations have been generated either at the shock or at the material interface. 

A considerably more difficult test that does not have an exact solution was suggested by [1^. This has a "slab" 
of helium/air mixture initially at rest in air. A moderate shock then strikes the slab, with the resulting shocks and 
rarefaction interacting multiple times. This is a Id restriction of the shock-bubble interaction used in the original 
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Ghost Fluid paper [TT!] . The relativistic analogue used here has initial data 




0.17933, 

0, 

0, 

0, 



1.57, 

P= 1, 
P= 1, 
1, 



7 = 1.4 

7 = 1.4 

7 = 1.67 

7 = 1.4 



X < 0.25 
0.25 < a; < 0.45 
0.45 <x < 0.55 

a; > 0.55. 



(4.2) 



The test is run to t = 0. 
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FIG. 3: The results for an isolated shock wave hitting a low density slab of material. Multiple reflected and transmitted waves 
can be seen, as expected. The movement and compression of the low density slab is also apparent. Velocity and pressure once 
again remain continuous at the interface as desired. TVD reconstruction using the MC limiter is used, and only 100 points are 
shown in each plot. 
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Results for this more complex interaction are shown in figure [3l TVD reconstruction using the MC limiter is used. 
The strong interacting features are captured, with the interfaces moving independently in reaction to the shocks. No 
oscillations occur, and the result converges with resolution. 



2. Perturbed shock 



A more complex example is that of a shock perturbed by a strong, smooth feature. This test, suggested by [26[ 
and used e.g. by [l^, does not have an exact solution, so is compared to a reference solution computed with 12800 
points. We modify this test to include an interface, using the initial data 

/ Po = 5, t; = 0, p = 50, 7 = 1.4 x< 0.5 , . 

\ po = 2 + 0.3sin(50a;), v = 0, p = 5, 7= 1.67 a; > 0.5 ' ^ ^ 

Again the domain is x e [0, 1] and the test is run to time t = 0.35. 

Results for this test are shown in figures d] and O It is clear that the results are converging to the correct solution 
irrespective of the numerical method used. No spurious oscillations are introduced. However, it appears that there 
is a slight discrepancy at the material interface. In particular, the edges of the smooth feature do not appear quite 
correct. We investigate this with another test. 



3. Moving sine wave 



A final flat space test is adapted from a Newtonian test suggested by [28| . A sine wave is advected in a surrounding 
fiow. The material interfaces should be stable, and so the evolution should be trivial. The initial data used here is 

Po= 1, v= 0.5, p= I, 7= 1.4 X < 0.16 

po= 1 + 0.3 sin(50(a;- 0.16)), 0.5, p= 1, 7= 1.67 0.16 < x < 0.537 . (4.4) 

po = 1, v= 0.5, p= 1, 7= 1.4 a; > 0.537 

Again the domain is a; € [0, 1] and the test is run to time t — 0.4. 

Results for this test are shown in figures [H] and [71 The majority of the features are well captured and no spurious 
oscillations are introduced at the material interfaces. It is clear, however, that the accuracy suffers near the interfaces 
due to the strong gradients there, and the overly simple Ghost Fluid boundary condition. This converges with 
resolution. We will see in the next section that in simple GR simulations modelling neutron stars this is not a major 
concern. 



B. GR tests 



In this section we look at tests in nonlinear general relativity, focusing on tests relevant for neutron stars. We use 
the system of equations outlined in section IIII A 21 All of the tests look at static or nonlinearly perturbed TOV stars 
in spherical symmetry. The initial data is always constructed with a polytropic EOS 

P = A>2 (4.5) 

where K is the polytropic constant, and evolved with the ideal fiuid EOS of equation (j3.5[) . where 7 may vary in 
space. 



1. Static stars 



As a reference solution we construct and evolve a standard static single component star. The initial data is 
determined by the initial central density po{t = 0, r = 0) denoted pc, and the EOS. Here we use the values 

Pc = 1.28xl0~^ 7 = 2, K=100. (4.6) 

We then evolve this using the TVD method as described in section IIII Bl The errors in the evolution are shown in 
figure [51 Second order convergence is seen over the majority of space at late time. However, near the surface of the 
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FIG. 4: Results for a test with a shock tube test with a sine wave density profile for 200 points. The higher accuracy when 
using PPM is evident in the post shock density profile, with the sine wave being captured more accurately. The solid line shows 
a reference solution computed using 12800 points. 

star convergence is lost, both due to the "clipping" inherent in HRSC methods, and due to the atmosphere algorithm 
used here. This is the dominant factor in reducing the order of convergence of the norms. 

As a first test of the accuracy of the Ghost Fluid algorithm in GR we then introduce a "trivial" interface into the 
TOV star at coordinate location r = 3.015. On either side of this interface the EOS is the same as that given in 
equation (|4.6p . We then evolve, imposing the Ghost Fluid boundary condition at the interface. It should be noted 
that in this simulation, although the level set is evolved the interface inferred from the zero of the level set remains 
within the same grid cell at all times. As can be seen from figure [9l by comparing to the reference solution of figure [51 
the Ghost Fluid method makes virtually no difference to the accuracy of the results. The convergence properties are 
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FIG. 5: Results for a test with a shock tube test with a sine wave density profile for 800 points. In addition to greater accuracy 
for the maxima and minima of the sine wave, similar improvements in accuracy at the shock wave are visible for this test. 



nearly identical, and the effect of the interface is very difficult to detect in the errors. 

To look at the effect of the Ghost Fluid boundary condition in more detail, in figure [TUl we show the difference in 
the radial 3-velocity component between the reference evolution without the interface and the evolution containing 
the interface. The effect of the interface is not directly obvious in the plot. The velocity should be continuous at 
the interface, and clearly is. The differences between the two evolutions appear near the surface of the star, which is 
where the velocity is largest. This indicates that the differences are due to minor reflections from the interface which 
are clearly at a very low level. In addition we look at the change in the entropy during the evolution. The entropy is 
extrapolated by the Ghost Fluid algorithm, so should show the effect of the boundary condition. The initial entropy is 
0(1), so whilst the effect of the boundary condition is clear in inducing a small jump at the interface, the magnitude 
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FIG. 6: An advected sine wave, evolved with a resolution of 200 points. Since the initial density profile is simply advected with 
the fluid velocity, an exact solution exists for this test. The differences between the two reconstruction methods are evident, 
with the PPM test has more accurately capturing the sine wave. 




FIG. 7: An advected sine wave, evolved with a resolution of 800 points. The sine wave is now captured well with both techniques, 
though at the interface between the fluids there is still some error. Comparing to figure |5] shows the (slow) convergence of this 
error. 



of the effect at ^ 10~ is clearly negligible. 

We tfien move on to look at static stars containing a genuine interface between two different components. We 
have studied a number of different configurations, with the qualitative behaviour being similar in all cases. As an 
illustration we use the initial data 

( Pc = 6x 10-4, 7=2, K= 100, r < 3.015 , . 

\ 7= |, X= 11.17, r> 3.015 " ^ ' 

This ensures the TOV is stable, and that pressure and its first derivative are continuous at the interface. This initial 
data is illustrated in figure [11] 

The result of the evolution is shown in figure [T^ In this case we do see an effect in the Hamiltonian constraint 
from the interface at r = 3.015. This leads to larger errors within the star (at fixed resolution). However, these errors 
appear to converge at second order (except for the points directly next to the interface). This means that, for these 
resolutions, the convergence of the constraint error is actually better than for the reference solution shown in figure [S] 
This is because the error at the surface, which does not converge at second order due to the atmosphere algorithm 
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FIG. 8: A single component TOV star is evolved as a reference case. The error in the Hamiltonian constraint is shown as a 
function of space at late time in the left panel, with the time-dependent behaviour of the 1-norm shown in the inset. Simulations 
with 640 and 1280 points are shown, with the behaviour in space scaled for second order convergence. The behaviour of the 
norm is scaled for convergence order 1.55. Second order convergence is observed nearly everywhere. The surface, located 
near r ~ 10, is the clear location where the order of convergence is reduced, as expected. This leads to the reduced order of 
convergence in the 1-norm. In the right panel the 1-norm of the density errors is shown, scaled for second order convergence. 
Here the impact of the surface errors is not so pronounced, so second order convergence in the norm is approached. 



Error scaled to second order with resolution, time = 1 000 Error scaled to second order with resolution 




Radius Time 

FIG. 9: A single component TOV star with an interface artificially introduced at r = 3.015 is evolved. These results should be 
compared with the reference solution in figure |8l The accuracy and convergence properties are nearly identical to the reference 
solution. 

and the inherent properties of HRSC methods, contributes less (relatively) to the total error. In contrast, the errors 
in the density are (relatively) larger at the surface as we are using a lower central density and a softer EOS for the 
outer fluid. Hence at late times as the surface errors increase the effects are seen in the errors of the density, which 
drifts away from perfect second order convergence. 
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FIG. 10: A closer look at the evolution of the single component TOV star with the artificial interface at r ~ 3. Only the 
simulation using 1280 points is shown. In the left panel the radial 3-velocity component is shown. As expected, the velocity is 
continuous across the interface. The difference between this evolution and the reference solution are due to interactions between 
the errors introduced by the surface and the interface, and are strongest near the surface. In the right panel the change in the 
entropy between the initial data and t ~ 270 is shown. As the entropy is extrapolated across the interface this shows the jump 
induced by the Ghost Fluid boundary condition; however, the size of this error is negligible. 
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FIG. 11: The initial data used for the evolution of a two component star. The outer fluid uses a softer EOS than the reference 
solution, leading to a larger radius. The interface at r ~ 3 is clearly seen from the density jump. 



2. Perturbed star 



The tests in section HV B II show that there is no problem in extending the level set method augmented with the 
Ghost Fluid boundary condition to GR when studying simple solutions. However, we really require that these methods 
will work in strongly nonlinear situations which will include interface motion and shock formation and interaction. In 
order to look at this in a 1+1 context we take a two component TOV star and give a strong nonlinear perturbation 
designed both to move the interface and to produce shocks propagating within the fluid. Due to the limitations of 
our atmosphere algorithm the results will not be trustworthy when the shocks start to reach the surface of the star, 
so we keep the evolutions relatively short. However, these evolutions are sufficient to validate the approach that we 
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FIG. 12: Errors introduced in evolving a two component star. This should be compared to the reference solution in figure |8] 
and the artificial interface solution in figure |9l The interface at r ~ 3 now has a visible effect in the errors, but these are 
still manageable and converge. Note that at these resolutions the convergence of the norm of the constraints is actually better 
than in the reference solution, as the errors in the interior (which converge at second order) are larger due to the non-trivial 
interface, whilst those near the surface (which converge at less than second order) are smaller due to the softer EOS. The errors 
in the density show good convergence until late times when the errors near the surface become sufficiently large. 



use here. 

The initial data used is based on a two component TOV star with parameters 

Pc = 1.28 X 10-3, 7=2, K= 100, r < 3.015 , . 

j= 1.9, K= 51.57, r> 3.015 " ^ 

The interior is then perturbed using 

Po = (/Oo)tov (1 + '^W) , P ^ Ptov {I + h{r)) , (4.9) 

where 

hir) = l^o(^-^^^'-m--m r<2.5, ^^^^^^ 
[0 r > 2.5. 

This is illustrated in the left panel of figure [131 The strong perturbation leads to a very steep feature near r ~ 2 that 
will form a shock and move the interface. As seen in the right panel of figure [121 there is an additional shock formed 
by the "ingoing" piece of the initial perturbation reflecting from the origin; at the time of this plot it has not reached 
the interface. 

The constraint errors in the evolution of the nonlinearly perturbed star are shown in figure 1141 As expected we 
now see a large number of features indicated on the plot. As in the reference solution in figure [SI we see the error 
from the surface, which is sizeable and does not converge cleanly at second order. As in the static two component 
test in figure [11] we see an error at the interface, which in this case is relatively small and converges except at a few 
points directly neighbouring the interface. Finally we see strong features at the two shocks where, due to the nature 
of HRSC methods, the order of accuracy is degraded to first order. 

Figure [T5l shows the motion of the interface and also shows that no spurious oscillations occur due to the interaction 
between the outgoing shock and the interface. The left panel shows the velocity at a time when the outgoing shock 
has passed through the interface located at r 3. The velocity is continuous with no spurious oscillations, despite 
the interaction with the shock. The right panel shows the motion of the interface. Initially it falls in, due to the 
additional matter in the interior. It then interacts with the shock, being pushed outwards and jumping between cells, 
before starting to fall back again, with the additional interaction with the reflected shock at t ~ 18 being insufficient 
to stop its inward motion. 
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Initial Data Time = 15 




Radius Radius 

FIG. 13: The left panel shows the initial data for the nonlinearly perturbed two component star. The interior perturbation 
near r ~ 2 is not quite a shock but will trigger an outgoing wave that steepens into a shock before passing through the interface 
at r ~ 3. It will also produce a strong ingoing wave that will reflect from the origin, producing another shock. The right panel 
shows the evolved at a time when the initial outgoing shock has passed through the interface but before the reflected shock has 
had time to reach it. 
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FIG. 14: Constraint errors in the evolution of the nonlinearly perturbed star. The curves are scaled for second order convergence, 
which is seen away from the surface and the shocks. 



18 




FIG. 15: The left panel shows the radial component of the 3- velocity for the nonlinearly perturbed star, indicating the locations 
of the shocks. By this point the interface, initially at r = 3.015, has moved outwards and has started to move back in again. 
As expected the velocity is smooth at the interface and no spurious oscillations have been generated by the interaction with 
the outgoing shock. The right panel indicates the interface location as a function of time, showing its coordinate location as 
inferred from the level set, and the computational coordinate of the centre of the cell within which the interface is found. We 
see that even with this strong perturbation the motion is small but similar for all resolutions. 



V. CONCLUSIONS 



In this paper we have considered models of relativistic matter containing sharp interfaces and their particular 
importance and implementation for numerical simulations. We expect that simulations of realistic models of neutron 
stars, as required for gravitational wave template calculations, will contain features that will be too sharp to be resolved 
as a smooth feature using reasonable numerical resolution. It follows that models incorporating sharp interfaces and 
techniques for evolving them cleanly will be required. 

The use of level set techniques to capture and evolve the interfaces that might form and disappear through mergers 
is already well known in numerical relativity, particularly in the context of event horizon finders. Their use here is 
natural and simple, and should allow the techniques investigated here to easily be extended to more spatial dimensions. 

However, the description of the interface location is not sufficient. We must impose a boundary condition at the 
interface, which will complete the description of the continuum model. The Ghost Fluid boundary condition used 
here has the advantage of simplicity, which makes it ideal for a proof-of-principle test such as this. However, it is 
questionable if this condition is sufficient for all situations, and more advanced techniques may be required with more 
complex EOSs than those considered here. 

The results of section HVl illustrate the advantages of the combined level set and Ghost Fluid method. It is simple, 
and for strong shocks or smooth flow interacting with interfaces it gives good results. We have shown that it directly 
extends to GR, and that for mildly perturbed neutron star tests the results are stable, convergent, and as accurate 
as the numerical method employed. However, when there is a long-lasting, strong feature that interacts with the 
interface, the test shown in section HV A 31 illustrates the limitations of the technique. This may indicate the need for 
better boundary conditions at the interface before high accuracy simulations for gravitational wave extraction can be 
attempted for complex neutron star models. 
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FIG. 16: A stable contact discontinuity, initially at a; = 0.5, with the fluid to the left having 7 = 5/3 and the fluid to the 
right having 7 = 4/3 is advected to the right with speed 0.1. The initial data is set such that the interface should propagate 
unchanged. The plot shows this initial data evolved using two resolutions and standard HRSC methods to time t = 0.2 using the 
colour function approach, with the exact solution given by the solid line. Spurious oscillations are introduced at the interface. 
These oscillations do not converge with resolution, indicating a failure of the approach. 



APPENDIX A: COLOUR FUNCTIONS 



At first glance it would appear that the simple tests shown in this paper are solved with an excessively complex 
numerical method. From the simplest point of view, our system of equations describing relativistic hydrodynamics 
has simply been augmented with an equation of motion for the EOS parameter 7. This can be described by a "colour 
function" and advected using standard HRSC methods if written in the form 

dt (pol) + {poiv') = 0. (Al) 

In fact any parameter n required by the EOS can be advected using an equation of this form. 

However, it is well known in the CFD literature that this naive approach fails for the simplest interfaces. Figure [TBI 
shows this approach failing for a stable material interface. The interface is set at a; = 0.5, with initial data 

/(1,0.1,|,§) x<0.5, 

leading to a stable interface propagating to the right. The standard conservation law form, as in section IIII A H 
augmented with the colour function equation (|A1[) . is evolved using the simple TVD HRSC method used in the rest 
of the paper. This is considerably simpler than the level set method, but as can be seen from the figure fails utterly 
to deal with this simple interface. Spurious oscillations are introduced. These oscillations are sizeable and do not 
converge with resolution. Clearly there is a failure of the approach, which is due to the formulation of the problem. 
The level set method combined with the Ghost Fluid boundary condition deals with this situation with no problem 
whatsoever, as illustrated in figure [T7l 

The reasons for the failure of the colour function approach have been studied in detail by a number of authors. 
As explained by [2§] the key is that the inherent "smearing" of the HRSC scheme will effectively "mix" the fluids 
either side of the interface in a narrow region. The detailed way that the mixing occurs will depend on the resolution 
and the numerical method, but it will always occur. However, without imposing further conditions it is clear that 
this mixing will not be thermodynamically consistent. That is, the numerical mixing will lead to spurious generation 
of entropy at the interface, creating the oscillations observed. Whilst it is possible to construct numerical methods 
that are thermodynamically consistent ( e.g . [1^; for example, for the EOS used here evolving the EOS parameter 
K = 1/(7 — 1) is consistent, as shown by \29\) this may not always be straightforward, and the problem of the length 
over which the interface is smeared, as described in the introduction, remains. 
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FIG. 17: The stable contact discontinuity test evolved to time t = 2 using the level set method with the Ghost Fluid boundary 
condition. Comparing to figure [16] we see that there are no spurious oscillations and the interface is captured perfectly at all 
resolutions. 
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